python_examples

Polyfem Tutorial

This is a jupyter notebook. The "real" notebook can be found here.

Polyfem relies 3 main objects:

  1. Settings that contains the main settings such discretization order (e.g., $P_1$ or $P_2$), material parameters, formulation, etc.
  2. Problem that describe the problem you want to solve, that is the boundary conditions and right-hand side. There are some predefined problems, such as DrivenCavity, or generic problems, such as GenericTensor.
  3. Solver that is the actual FEM solver.

The usage of specific problems is indented for benchmarking, in general you want to use the GenericTensor for tensor-based PDEs (e.g., elasticity) or GenericScalar for scalar PDEs (e.g., Poisson).

A typical use of Polyfem is:

settings = polyfempy.Settings()
# set necessary settings
# e.g. settings.discr_order = 2

problem = polyfempy.GenericTensor() # or any other problem
# set problem related data
# e.g. problem.set_displacement(1, [0, 0], [True, False])

settings.set_problem(problem)

#now we can create a solver and solve
solver = polyfempy.Solver()

solver.settings(settings)
solver.load_mesh_from_path(mesh_path)

solver.solve()

Note 1: for lecacy reasons Polyfem always normalizes the mesh (i.e., rescale it to lay in the $[0,1]^d$ box, you can use setting.normalize_mesh = False to disable this feature.

Note 2: the solution $u(x)$ of a FEM solver are the coefficients $u_i$ you need to multiply the bases $\varphi_i(x)$ with: $$ u(x)=\sum u_i \varphi_i(x). $$ The coefficients $u_i$ are unrelated with the mesh vertices because of reordering of the nodes or high-order bases. For instance $P_2$ bases have additional nodes on the edges which do not exist in the mesh.

For this reason Polyfem uses a visualization mesh where the solution is sampled at the vertices. This mesh has two advantages:

  1. it solves the problem of nodes reordering and additional nodes in the same way
  2. it provides a "true" visualization for high order solution by densely sampling each element (a $P_2$ solution is a piecewise quadratic function which is visualized in a picewise linear fashion, thus the need of a dense element sampling).

To control the resolution of the visualization mesh use settings.vismesh_rel_area.

Examples

Some imports for plotting

In [1]:
import plotly.offline as plotly
import plotly.graph_objs as go
import plotly.figure_factory as ff

#Necessary for the notebook
plotly.init_notebook_mode(connected=True)

algebra

In [2]:
import numpy as np

stuff for the animation

In [3]:
import ipywidgets as widgets
from ipywidgets import interact

and finallypolyfempy

In [4]:
import polyfempy as pf

Plotting utility

This code is not particularly interesting.

It converts a tet-mesh (4 indices) into a face-based tri mesh and uses plotly Mesh3d to plot it.

In [5]:
def create_plot_mesh_and_layout(vertices, connectivity, function):
    x = vertices[:,0]
    y = vertices[:,1]

    if vertices.shape[1] == 3:
        z = vertices[:,2]
    else:
        z = np.zeros(x.shape, dtype=x.dtype)

    if connectivity.shape[1] == 3:
        f = connectivity
    else:
        #Convert a tet-mesh into face based triangles.
        f = np.ndarray([len(connectivity)*4, 3], dtype=np.int64)
        for i in range(len(connectivity)):
            f[i*4+0] = np.array([connectivity[i][1], connectivity[i][0], connectivity[i][2]])
            f[i*4+1] = np.array([connectivity[i][0], connectivity[i][1], connectivity[i][3]])
            f[i*4+2] = np.array([connectivity[i][1], connectivity[i][2], connectivity[i][3]])
            f[i*4+3] = np.array([connectivity[i][2], connectivity[i][0], connectivity[i][3]])

    mesh = go.Mesh3d(x=x, y=y, z=z,
                     i=f[:,0], j=f[:,1], k=f[:,2],
                     intensity=function, flatshading=function is not None,
                     colorscale='Viridis',
                     contour=dict(show=True, color='#fff'),
                     hoverinfo="all")
    layout = go.Layout(scene=go.layout.Scene(
        aspectmode='data',

        xaxis=dict(
            autorange=True,
            showgrid=False,
            zeroline=False,
            showline=False,
            ticks='',
            showticklabels=False
        ),
        yaxis=dict(
            autorange=True,
            showgrid=False,
            zeroline=False,
            showline=False,
            ticks='',
            showticklabels=False
        ),
        zaxis=dict(
            autorange=True,
            showgrid=False,
            zeroline=False,
            showline=False,
            ticks='',
            showticklabels=False
        )
    ))

    return mesh, layout
In [6]:
def plot(vertices, connectivity, function, camera=None):
    mesh, layout = create_plot_mesh_and_layout(vertices, connectivity, function)
    if camera is not None:
        layout.scene.camera = camera

    fig = go.Figure(data=[mesh], layout=layout)

    plotly.iplot(fig)

Creates a quad mesh of n_pts x n_pts in the form of a regular grid

In [7]:
def create_quad_mesh(n_pts):
    extend = np.linspace(0,1,n_pts)
    x, y = np.meshgrid(extend, extend, sparse=False, indexing='xy')
    pts = np.column_stack((x.ravel(), y.ravel()))

    faces = np.ndarray([(n_pts-1)**2, 4],dtype=np.int32)

    index = 0
    for i in range(n_pts-1):
        for j in range(n_pts-1):
            faces[index, :] = np.array([
                j + i * n_pts,
                j+1 + i * n_pts,
                j+1 + (i+1) * n_pts,
                j + (i+1) * n_pts
            ])
            index = index + 1

    return pts, faces

Plate hole

This is the python version of the plate with hole example explained here.

Set the mesh path

In [8]:
mesh_path = "plane_hole.obj"

create settings

In [9]:
settings = pf.Settings()

pick linear $P_1$ elements. If the mesh would be a quad it would be $Q_1$

In [10]:
settings.discr_order = 1

normalize the mesh to be in $[0,1]^2$

In [11]:
settings.normalize_mesh = True

and choose Young's modulus and poisson ratio

In [12]:
settings.set_material_params("E", 210000)
settings.set_material_params("nu", 0.3)

We are use a linear material model

In [13]:
settings.tensor_formulation = pf.TensorFormulations.LinearElasticity

Next we setup the problem

In [14]:
problem = pf.GenericTensor()

sideset 1 has zero displacement in $x$

In [15]:
problem.set_displacement(1, [0, 0], [True, False])

sideset 4 has zero displacement in $y$

In [16]:
problem.set_displacement(4, [0, 0], [False, True])

sideset 3 has a force of [100, 0] applied

In [17]:
problem.set_force(3, [100, 0])

fianally set the problem

In [18]:
settings.set_problem(problem)

Solve!

In [19]:
solver = pf.Solver()

solver.settings(settings)
solver.load_mesh_from_path(mesh_path)

solver.solve()

Get the solution

In [20]:
[pts, tets, disp] = solver.get_sampled_solution()

diplace the mesh

In [21]:
vertices = pts + disp

and get the stresses

In [22]:
mises, _ = solver.get_sampled_mises_avg()

finally plot with the above code

In [23]:
top_camera = dict(
    up=dict(x=0, y=1, z=0),
    center=dict(x=0, y=0, z=0),
    eye=dict(x=0, y=0, z=0.8)
)

plot(vertices, tets, mises, top_camera)

Note that we used get_sampled_mises_avg to get the Von Mises stresses because they depend on the Jacobian of the displacement which, becase we use $P_1$ elements, is piece-wise constant. To avoid that effect in get_sampled_mises_avg the mises are averaged per vertex weighted by the area of the triangles. If you want the "real" mises just call

In [24]:
mises = solver.get_sampled_mises()
plot(vertices, tets, mises, top_camera)

Torsion

Non-linear example. We want to torque a 3D bar around the $z$ direction.

The example is really similar as the one just above.

Sets the mesh and create the settings. As before

In [25]:
mesh_path = "square_beam_h.HYBRID"

settings = pf.Settings()

It is an hex-mesh so we are using $Q_1$

In [26]:
settings.discr_order = 1

In this case the mesh has already the correct size.

In [27]:
settings.normalize_mesh = False

Choose Young's modulus and poisson ratio, as before

In [28]:
settings.set_material_params("E", 200)
settings.set_material_params("nu", 0.35)

Differently from before we want a non-linear material model: NeoHookean

In [29]:
settings.tensor_formulation = pf.TensorFormulations.NeoHookean

and we want to do 5 steps of incremental loading to avoid ambiguities in the rotation

In [30]:
settings.nl_solver_rhs_steps = 5

Now we setup problem with fixed sideset, rotating an number of tours

In [31]:
problem = pf.Torsion()

sideset 5 is fixed

In [32]:
problem.fixed_boundary = 5

sideset 6 rotates

In [33]:
problem.turning_boundary = 6

around the $z$-axis (2)

In [34]:
problem.axis_coordiante = 2

by half a tour

In [35]:
problem.n_turns = 0.5

Now we choose a coarse visualization mesh

In [36]:
settings.vismesh_rel_area = 0.001

and set the problem and solve

In [37]:
settings.set_problem(problem)

#solving!
solver = pf.Solver()

solver.settings(settings)
solver.load_mesh_from_path(mesh_path)

solver.solve()

takes approx 1 min because it is a complicated non-linear problem!

Get solution and stesses as before

Since we want to show only the surface there is no need of getting the whole volume, so we set boundary_only to True

In [38]:
[pts, tets, disp] = solver.get_sampled_solution(boundary_only=True)
vertices = pts + disp
mises, _ = solver.get_sampled_mises_avg(boundary_only=True)

and plot the 3d result!

In [39]:
plot(vertices, tets, mises)

Fluid simulation

Create the mesh using the utility function

In [40]:
pts, faces = create_quad_mesh(50)

create settings

In [41]:
settings = pf.Settings()

settings.vismesh_rel_area = 0.001

pick linear $Q_2$ elements for velocity and $Q_1$ for pressure

In [42]:
settings.discr_order = 2
settings.pressure_discr_order = 1

Set the viscosity of the fluid

In [43]:
settings.set_material_params("viscosity", 1)

We select stokes as material model

In [44]:
settings.tensor_formulation = pf.TensorFormulations.Stokes

The default solver do not support mixed formulation, we need to choose UmfPackLU

In [45]:
settings.set_advanced_option("solver_type", "Eigen::UmfPackLU")

We use the standard Driven Cavity problem

In [46]:
problem = pf.DrivenCavity()

we set the problem

In [47]:
settings.set_problem(problem)

We create the solver and set the settings

In [48]:
solver = pf.Solver()
solver.settings(settings)

This time we are using pts and faces instead of loading from the disk

In [49]:
solver.set_mesh(pts, faces)

Solve!

In [50]:
solver.solve()

We now get the solution and the pressure

In [51]:
[pts, tris, velocity] = solver.get_sampled_solution()

and plot it!

In [52]:
top_camera = dict(
    up=dict(x=0, y=1, z=0),
    center=dict(x=0, y=0, z=0),
    eye=dict(x=0, y=0, z=1.2)
)

plot(pts, tris, np.linalg.norm(velocity, axis=1), top_camera)
In [53]:
n_pts = len(pts)
scaling = 3
quiver = ff.create_quiver(
    pts[0:n_pts:10,0], pts[0:n_pts:10,1],
    scaling*velocity[0:n_pts:10,0], scaling*velocity[0:n_pts:10,1])

quiver.layout.yaxis.scaleanchor="x"
quiver.layout.yaxis.scaleratio=1
plotly.iplot(quiver)

Time dependent simulation

Create the mesh using the utility function

In [54]:
pts, faces = create_quad_mesh(50)

create settings

In [55]:
settings = pf.Settings()

pick linear $Q_1$ elements.

In [56]:
settings.discr_order = 1

and choose Young's modulus and poisson ratio

In [57]:
settings.set_material_params("E", 1)
settings.set_material_params("nu", 0.3)

We are use a linear material model

In [58]:
settings.tensor_formulation = pf.TensorFormulations.LinearElasticity

For efficienty in the browser we select a coarse vis mesh

In [59]:
settings.vismesh_rel_area = 0.001

We simulate from 0 to 10s and 50 steps.

In [60]:
settings.tend = 10
settings.time_steps = 50

Next we setup the problem, this doesnt have any parameters. It will...

In [61]:
problem = pf.Gravity()

we set the problem

In [62]:
settings.set_problem(problem)

We create the solver and set the settings

In [63]:
solver = pf.Solver()
solver.settings(settings)

This time we are using pts and faces instead of loading from the disk

In [64]:
solver.set_mesh(pts, faces)

Solve!

In [65]:
solver.solve()

Get the solution and the mises

In [66]:
pts = solver.get_sampled_points_frames()
tris = solver.get_sampled_connectivity_frames()
disp = solver.get_sampled_solution_frames()
mises = solver.get_sampled_mises_avg_frames()

Now we are ready to do the animation

First create a figure widget

In [67]:
top_camera = dict(
    up=dict(x=1, y=0, z=0),
    center=dict(x=0, y=0, z=0),
    eye=dict(x=0, y=0, z=1.2)
)

mesh, layout = create_plot_mesh_and_layout(pts[-1],tris[-1],mises[-1])
mesh.cmin = 0
mesh.cmax = 0.25
layout = go.Layout(scene=go.layout.Scene(
        camera=top_camera,

        aspectratio = dict( x=1, y=1, z=1 ),
        aspectmode = 'manual',

        xaxis = dict(range=[-0.1, 1.1]),
        yaxis = dict(range=[-0.1, 1])
))
fig = go.FigureWidget(data=[mesh], layout=layout)

then we need the callback

In [68]:
def on_value_change(value):
    frame_index = value['new']

    fig.data[0].x=pts[frame_index][:,0] + disp[frame_index][:,0]
    fig.data[0].y=pts[frame_index][:,1] + disp[frame_index][:,1]

    fig.data[0].intensity = mises[frame_index]

finally we create an animation widged (works only in the notebook)

In [69]:
play = widgets.Play(
value=0,
min=0,
max=len(mises)-1,
step=1,
description="Press play",
disabled=False
)

slider = widgets.IntSlider(min=0, max=len(mises)-1)

slider.observe(on_value_change, 'value')
play.observe(play, 'value')

widgets.jslink((play, 'value'), (slider, 'value'))
widgets.VBox((widgets.HBox((play, slider)), fig))